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MOLECULAR  DYNAMICS  SIMULATION 
OF  SEDIMENTATION 

R.A.J.  Borg,  D.A.  Jones  and  S.G.  Lambrakos* 

Materials  Research  Laboratory,  D.S.T.O.,  Melbourne,  Australia  SOSt 
*  Naval  Research  Laboratory,  Washington,  D.C.,  USA 

1.  INTRODUCTION 

Composition  B  is  a  widely  used  military  explosive  consisting  of  a  60/40  mixture  of  RDX 
(hexahydrotrinitro-s-triazine)  and  TNT  (trinitrotoluene).  Ordnance  is  filled  by  a  casting 
process  in  which  Composition  B  is  heated  above  the  melting  point  of  TNT  (82®C),  poured 
into  the  shell,  and  allowed  to  cool  and  solidify.  During  this  process  the  heated  mixture  is 
essentially  a  suspension  of  solid  RDX  particles  in  liquid  TNT,  and  sedimentation  of  RDX 
can  occur  in  the  shell  during  cooling  and  solidification  of  the  TNT.  The  sedimentation  will 
lead  to  localized  regions  of  higher  RDX  content  and  higher  sensitivity,  because  RDX  is  more 
sensitive  to  initiation  than  TNT. 

This  paper  describes  the  application  of  a  Molecular  Dynamics  (MD)  computer  simulation 
technique  which  was  used  to  study  the  sedimentation  process.  We  considered  N  macroscopic 
size  spherical  particles  in  a  viscous  fluid  subject  to  gravity  and  interparticle  forces.  Calcu¬ 
lations  were  performed  for  three  different  particle  sizes  for  both  monomodal  and  bimodal 
particle  size  distributions.  The  monomodal  results  showed  a  size  dependent  packing  effect, 
while  with  bimodal  distributions  fractionation  was  found  to  occur,  giving  rise  to  a  layer  of 
smaller  particles  on  the  top  of  the  sediment. 

The  CPU  time  for  the  initial  version  of  our  code  varied  as  because  of  the  nearest 
neighbours  problem.  This  limited  the  maximum  number  of  particles  in  our  simulations  to 
slightly  more  than  200,  and  so  later  versions  of  the  c.  Je  included  a  Monotonic  Lagrangian  Grid 
(MLG)  algorithm  (Boris,  1986)  to  keep  track  of  nearest  neighbour  relationships  and  reduce 
the  computational  time  to  an  order  N  dependence.  We  briefly  describe  the  implementation 
of  the  MLG  into  our  sedimentation  code  and  the  improvement  in  performance  which  was 
obtained. 

2.  MOLECULAR  DYNAMICS  MODEL 

In  MD  calculations  the  forces  on  all  particles  are  calculated  and  the  particles  we  moved 
according  to  Newton’s  equations  of  motion.  The  initial  state  of  the  sedimenting  system 
consists  of  N  particles  dispersed  homogeneously  in  a  viscous  medium  and  subject  to  gravity 
and  interparticle  forces.  The  final  state  consists  of  a  completely  sedimented  system  in  which 
none  of  the  particles  has  any  significant  motion. 

The  important  forces  in  a  sedimentation  process  ve  viscosity,  gravity,  interparticle  forces 
and  particle/floor  forces.  The  magnitude  of  the  viscous  force  on  a  spherical  particle  is  given 
by 

=  -67rr,riVi  (1) 

where  rj  is  the  liquid  phase  viscosity,  is  the  radius  of  particle  i  and  Vj  the  velocity.  This 
force  acts  in  a  direction  opposite  to  the  direction  of  motion  of  the  particle. 


•  •  • 
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€  X  10*^ 

Particle  Radi 

us  (/im) 

{J) 

30 

50 

100 

30 

0.473 

1.50 

10.4 

50 

1.50 

3.65 

18.5 

100 

10.4 

18.5 

58.4 

Table  1:  f  values  used  in  the  interparticle  force  calculation. 


<T  X  10* 

Particle  Radius 

(m) 

30 

50 

100 

30 

14.7 

19.6 

31.9 

50 

19.6 

24.5 

36.8 

100 

31.9 

36.8 

49.0 

Table  2:  tr  values  used  in  the  interparticle  force  calculation. 


Rad'us 

(pm) 

30 

50 

100 

K 

(V.m'*) 

3.54  X  10“” 

1.25  X  10-** 

8.22  X  10~*® 

Table  3:  K  values  used  in  the  particle/fioor  force  calculation. 


The  interparticle  forces  describe  the  interaction  between  particles.  Here,  the  particles 
are  assumed  to  be  spherical  "soft”  spheres.  A  soft  sphere  model  is  necessary  to  overcome 
numerical  stability  problems  as  well  as  reducing  particle  overlap.  A  truncated  Lennard-Jones 
12-6  type  potential  proved  to  be  the  most  suitable  and  the  force  is  given  by 


4ei 


12W/ 

,13 


(2) 


where  Vij  is  the  interparticle  distance,  tij  and  Vij  are  constants  for  pair  (ij)  and  d  is  a  unit 
vector  along  the  interparticle  line. 

A  particle-floor  force  is  needed  between  the  particles  and  the  bottom  of  the  container. 
The  form  of  this  force  was  taken  to  be 


JffioOT  _  Kj 

* «  ~  “ 
r  • 


(3) 


where  is  the  distance  between  particle  i  and  the  bottom  of  the  contuner,  Ki  is  a 

constant  for  particle  type  i  and  k  is  a  unit  vector  in  the  Z  direction.  The  form  of  this 
potential  is  basically  the  repulsive  part  of  Equation  (2)  and  was  chosen  for  consistency,  as 
well  as  its  success  in  providing  an  adequate  description  of  the  bottom  boundary. 

Values  of  eij  and  ffij  can  be  obtained  in  the  literature  for  atoms  and  molecules.  However, 
values  for  macroscopic  particles  like  those  used  in  this  study  are  not  readily  available.  The 
values  used  here  have  been  chosen  somewhat  arbitrarily  but  satisfy  the  required  conditions 
well,  i.e.  the  particles  behave  as  soft  spheres.  K  values  are  also  assigned  in  a  similar  fashion. 
Tables  1,  2  and  3  list  values  of  €ij,  Cij  and  Ki  used  in  this  work.  The  results  we  obtained 
were  found  not  to  be  sensitive  to  the  values  used. 
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The  effect  of  Brownian  motion  on  the  sedimenting  particles  has  not  been  included  as  the 
particle  sizes  considered  here  are  large  enough  for  it  to  be  neglected.  A  simple  criterion  for 
the  neglect  of  Brownian  motion  is  that  the  particle  Peclet  number  be  much  larger  than  unity 
(Bossis  and  Brady,  1984).  For  the  system  considered  here  this  will  be  the  case  provided  the 
particle  number  is  greater  than  1pm.  The  smallest  particle  radius  considered  here  is  30pm, 
and  so  Brownian  motion  has  not  been  included  in  the  calculation. 

The  total  force  on  a  particle  is  therefore  given  by 

N 

Fi  =  miBi  =  ^  Fij  +  rriig  -  6irjjriVi  +  (4) 


where  N  is  the  number  of  particles,  Vj  is  the  velocity  of  particle  i  and  the  acceleration. 
We  integrate  the  equations  of  motion  using  the  Verlet  equations  in  the  following  form 

Vi{0  =  ^  (ri(<  +  At)  -  ri(t  -  At))  (5) 

ri(t  +  At)  =  2ri(t)  -  ri(t  -  At)  +  a^(t)At*  (6) 

To  cidculate  Ti{t  +  At),  Bi{t)  must  be  known,  and  to  calculate  ai(t)  from  Equation  (4)  we 
need  Vi(t)  which,  in  turn,  is  dependent  on  ri(t  +  At).  This  difficulty  can  be  overcome  as 
follows;  Equation  (4)  can  be  written  (in  one  dimension)  as 


Ai-ai- 


•  • 


where 


1  /  ^  \ 

'  VciSo  / 


Equations  (7)  and  (8)  can  then  be  used  to  simplify  Equation  (5)  to  the  following  form 


*"  =  (.  +  T j  T' y  j 


Equation  (10)  can  now  be  evaluated  using  only  Xi*  and  Xj**  is  approximated  via  : 

*<‘=Xi»  +  ^At  (11) 

Each  time  step  integration  then  cycles  through  the  following  sequence:  evaluate  A  from 
Equation  (8),  solve  Equation  (10)  to  give  the  particle  velocities,  solve  Equation  (7)  to  give 
the  particle  accelerations,  then  solve  Equation  (6)  to  find  position  at  the  new  time  step.  In 
this  scheme,  the  force  evaluation  and  solution  of  the  equations  of  motion  are  interleaved  to 
overcome  the  velocity  dependence  of  the  force. 


•  • 
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Figure  1.  Top  surface  of  the  sediment  approximated  by  a  grid  of  squares. 


Initially,  the  particles  are  placed  randomly  within  a  rectangular  box  and  assigned  random 
velocities.  Periodic  boundary  conditions  are  applied  in  the  X  and  Y  directions,  and  gravity 
acts  in  the  negative  Z  direction.  The  number  of  particles  in  the  simulations  varied  from 
N  =  125  to  N=1000,  and  a  typical  valqe  for  the  time  step  was  200ps.  Complete  sedimentation 
usually  occured  within  approximately  12,000  to  14,000  cycles.  The  initial  version  of  the  code 
used  a  lookup  table  to  evaluate  Equation  (2),  but  the  MLG  version  calculated  the  force  at 
each  time  step. 

In  order  to  compare  results  for  different  particle  size  distributions  a  method  for  calculat¬ 
ing  the  sediment  volume  was  required.  Given  the  box  dimensions  and  the  particle  coordinates 
and  radii,  the  volume  can  be  calculated  by  evaluating  the  volume  under  the  top  surface  of 
the  particles,  and  this  can  be  defined  by  lowering  a  grid  of  small  squares  onto  the  sediment 
until  aU  of  the  squares  have  intersected  a  particle.  The  volume  is  then  the  sum  of  the  volumes 
under  every  square.  Figure  (1)  shows  a  typical  picture  of  the  top  surface  of  a  sediment  which 
has  been  approximated  using  this  method. 

The  imtial  version  of  the  program  was  designed  to  run  on  a  VAX  8700  computer  and 
practical  limitations  restricted  the  number  of  particles  in  the  simulation  to  around  N=200 
because  the  CPU  time  varied  as  N^.  This  is  a  common  problem  in  MD  simulations  and 
arises  because  in  a  system  of  N  particles  randomly  distributed  in  space  the  motion  of  any 
one  psirticle  is  in  principle  determined  by  the  remaining  N-1  other  particles.  In  practice 
however  the  particle  will  be  influenced  significantly  by  only  a  relatively  small  number  of 
nearby  particles,  defined  to  be  the  particles  ’’nearest  neighbours”,  and  finding  these  nearest 
neighbours  in  an  efficient  manner  is  a  central  part  of  any  MD  simulation. 

We  redesigned  our  code  to  overcome  this  problem  by  implementing  the  MLG  algorithm  to 
track  near  neighbour  relationships.  The  MLG  is  a  data  structure  in  which  adjacent  particles 
in  space  have  close  grid  indices.  The  data  structure  is  arranged  so  that  a  particle’s  nearest 
neighbours  are  readily  identified,  but  without  the  need  to  continually  calculate  the  distance 
to  each  of  the  other  N-1  particles  (an  order  N  operation).  The  computational  coat  of  the 
scheme  scales  as  N,  and  the  algorithm  is  ideally  suited  to  vectorization  because  it  uses  data 
&om  contiguous  memory  locations.  The  MLG  has  previously  been  used  for  a  variety  of  MD 
simulations  of  atomic  and  molecular  systems  (Lambrakos  et  al,  1989a,b),  but  this  is  the  first 
time  it  has  been  applied  to  macroscopic  size  particles  in  a  sedimenting  system. 


Stiimentatton 


131 


30/im 


50/xm 


HEIGHT  itim) 

Figure  3a.  Probability  of  finding  a  particle  as  a  function  of 
height  above  the  fioor  of  the  container  for  the  SO/50/tm  simulation. 

3.  RESULTS  AND  DISCUSSION 

Initial  calculationx  were  performed  for  suspensions  in  which  all  particles  had  the  same  radius. 
Three  particle  sises  were  investigated;  r^SO,  50  and  100pm,  for  N=125  or  216.  In  order  to  ob¬ 
tain  statistically  meaningful  results  six  replicate  calculations  were  performed  for  eadi  radius, 
the  only  difference  between  the  calculations  being  in  the  starting  positions  and  vdodties.  A 
relative  density  pr<  defined  as  the  ratio  of  sediment  density  to  particle  material  density,  was 
calculated  from  the  averaged  results  for  each  radius  value. 

If  the  sedimented  nwterial  had  sufficient  time  to  adjust  to  an  optimum  close  packed 
struettrre  we  would  find  pr  =  0.74,  independent  of  particle  use,  whereas  the  siinnlations  gave 
pr  =  0.72,  0.70,  0.68  for  particle  radius  values  of  30,  50,  and  100pm  respectively.  Thus  the 
packing  density  was  found  to  be  particle  sise  dependent,  with  an  increase  in  particle  sise 
giving  a  decrease  in  particle  density. 

One  possible  explanation  for  this  lies  in  the  terminal  velocities  of  the  particles.  By 
considering  the  free  fall  of  a  particle  in  the  suspension  (free  fall  here  means  that  there  are  no 
interparticle  or  particle/fioor  forces  on  the  particle)  it  can  be  shown  that  the  terminal  velocity 
of  a  particle  is  proportional  to  the  square  of  the  particle  radiiu.  Hence  the  smaller  particles 
have  a  slower  terminal  velocity  and  therefore  have  more  time  to  rearrange  and  adopt  a  more 
efficient  structure  than  the  larger  panicles. 

Three  different  bimodal  particle  sise  distributions  were  investigated;  30/50,  50/100,  and 
30/100pm  mixes.  With  bimodal  distributions  there  is  the  possibility  of  fractionation  occuring 
during  sedimentation,  and  therefore  plots  of  the  probability  of  locating  a  particle  centre  as  a 
function  of  height  (in  the  s  direction)  were  obtained.  These  are  shown  in  Figure  (2). 

A  cmnmon  feature  of  these  plots  is  a  pair  of  well  defined,  sharp  peaks  at  low  heights.  Fw 
example,  the  30/50pm  plot  (Figure  (2a))  has  a  sharp  peak  at  26fun  for  the  30pm  partides 
and  another  at  46pm  for  the  50pm  particles.  These  peaks  correspond  to  particles  in  contact 
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HEIGHT  {urn) 

Figure  3b.  Probability  of  finding  a  particle  as  a  function  of 
height  above  the  fioor  of  the  container  for  the  SO/lOOpm  siinnlation. 
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Figure  3c.  Probability  of  finding  a  particle  ae  a  ftinciion  of 
height  above  the  floor  of  the  container  for  the  SO/lOOpm  nnnilation. 


•  •  • 


•  • 
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log  N 

Figure  3.  CPU  time  for  100  cycles  versus 
total  number  of  particles  in  the  simulation. 


with  the  floor  of  the  container.  The  fact  that  these  peaks  occiir  at  heights  slightly  lower  than 
the  corresponding  particle  radius  indicate  a  small  amount  of  overlap  between  the  floor  and 
the  particles  (i.e.  squashing). 

The  tendency  for  the  smaller  particles  to  settle  after  the  larger  particles  is  shown  by  the 
plots.  For  the  30/50/im  simulation  the  30/im  probability  curve  extends  out  to  360/im  whereas 
the  SOpm  curve  only  goes  out  to  310pm.  Thus,  beyond  310pm  there  is  zero  probability  of 
finding  a  50pm  particle  and  only  30pm  particles  will  be  found  in  the  range  310pm  to  360pm, 
i.e.  a  layer  of  30pm  particles  will  be  formed  on  the  top  of  the  sediment.  Similarly,  for  the 
SO/lOOpm  simulation  (Figure  (2b))  100pm  particles  are  not  found  above  578pm  and  only 
50pm  particles  are  found  between  578pm  and  623pm.  The  30/100pm  simulation  (Figure 
(2c))  is  not  as  clear  cut  since  the  probability  of  finding  a  100pm  particle  above  470pm  is  not 
zero,  but  since  it  is  very  small  we  can  take  470pm  as  the  upper  limit  for  the  100pm  particles 
here.  Only  30pm  particles  will  be  found  between  470pm  and  530pm. 

The  bimodal  results  were  obtained  using  either  N=216  or  N=512,  with  the  larger  number 
of  particles  being  rim  on  the  MLG  version  of  the  code.  We  also  ran  a  series  of  simulations 
using  up  to  1000  particles  to  check  the  order  N  scaling  of  the  MLG  on  both  the  VAX  8700  and 
a  Cray  X-MP.  Plots  of  CPU  time  for  100  cycles  versus  total  number  of  particles  are  shown  in 
Figure  (3)  for  both  the  oripnal  version  of  the  code  and  the  MLG  version  on  both  the  VAX 
and  the  Cray. 

The  original  code  shows  the  dependency  (gradient  equals  1.98),  while  the  MLG 
code  on  the  VAX  shows  a  linear  dependence  on  N  (gradient  equals  1.08),  as  expected.  It  is 
interesting  to  note  that  the  MLG  code  on  the  VAX  takes  longer  to  run  100  cycles  than  does 
the  original  version  of  the  code  when  small  particle  numbers  are  used.  This  is  not  unexpected 
because  of  the  additional  initial  computational  overheads  associated  with  the  MLG,  and  the 
repeated  use  of  Equation  (2)  instead  of  a  lookup  table.  Crossover  occurs  around  N=500,  and 
for  N=:1000  the  MLG  version  is  more  than  twice  as  fast. 

The  MLG  code  when  run  on  the  Cray  also  shows  a  linear  dependence  on  N  (gradient 
equals  0.97)  and  considerably  reduced  CPU  times.  The  speedup  is  only  marginally  better 
than  that  due  to  the  faster  clock  rate  on  the  Cray  however,  indicating  that  little  vectorization 
is  occuring,  and  this  is  due  to  the  relatively  small  number  of  particles  we  are  using.  In  an 


•  • 
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MLG  code  the  N  particles  eure  specified  by  assigning  values  to  three  numbers:  NX,  NY,  and 
NZ,  where  N  =  NX  x  NY  x  NZ.  In  a  "regular”  MLG  structure  the  computationally  intensive 
parts  of  the  code  are  contained  within  a  nested  set  of  DO  loops  over  indices  i,  j,  and  k,  which 
vary  between  1  and  NX,  NY,  and  NZ  respectively.  For  1000  particles  we  have  NX  =  NY  = 
NZ  =10,  which  means  that  the  innermost  loop  has  a  count  of  only  10,  which  is  far  too  short 
for  the  increase  in  speed  due  to  the  pipeline  architecture  to  be  effective. 

4.  CONCLUSION 

For  monomodal  puticie  size  distributions  the  MD  results  show  that  the  packing  density  is 
particle  size  dependent,  with  an  increase  in  particle  radius  giving  a  decrease  in  sediment 
density.  For  bimodal  distributions  fractionation  was  found  to  occur  and  a  layer  of  smaller 
particles  formed  on  the  top  of  the  sediment.  This  type  of  inhomogeneity  has  important  con¬ 
sequences  for  the  sensitivity  and  mechanical  properties  of  an  explosive  where  sedimentation 
may  have  occurred  during  production.  Use  of  the  regular  MLG  algorithm  reduced  the  N^ 
dependency  on  CPU  time  to  order  N,  but  further  improvements  in  efficiency  require  the  use 
of  either  larger  particle  numbers  or  a  "skew- periodic”  MLG  on  a  vector  compiler  (Lambrakos 
and  Boris,  1987),  or  the  application  of  the  regular  MLG  on  a  massively  parallel  computer 
(Phillips,  1991). 
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